
#this needs to happen at the kecamatan, rather than village level, since that's the lowest administrative unit for which we have mosque data

kec = readShapePoly("./_3_data/map_data/admin_maps/gadm36_IDN_shp/gadm36_IDN_3.shp")
kec@data$id = rownames(kec@data)
kec.points = fortify(kec)
kec@data$prov_id = substring(kec@data$CC_3, first = 1, last = 2)
kec = kec[!is.na(kec@data$prov_id),]
kec = kec[kec@data$prov_id > 30 & kec@data$prov_id < 40,]

op = readShapePoly("./_3_data/map_data/opium_maps/opium_banned_1886.shp")
boundary = readShapeLines("./_3_data/map_data/opium_maps/opium_boundary_1886.shp") 

mosques_data = read_csv("./_3_data/mosque_data/mosques_mush_fin_9417.csv")

#intersection
clip = gIntersection(op, kec, byid = TRUE, drop_lower_td = TRUE) #clip polygon 2 with polygon 1
clip.points = fortify(clip)

#calculate kecamatan distance to boundary
df <- matrix(NA, ncol = 5, nrow = 30000)

for(i in 1:length(kec@polygons)){
  df[i,] <- c(i,
              as.character(dist2Line(kec@polygons[[i]]@Polygons[[1]]@labpt, boundary)[1,1]), 
              as.character(kec@polygons[[i]]@Polygons[[1]]@labpt),
              as.character(kec@polygons[[i]]@ID))
}



distances =
  df %>%
  data.frame() %>%
  set_colnames(c("rowid", "distance", "long", "lat", "id")) %>%
  dplyr::select(-rowid) %>%
  na.omit() %>%
  mutate_at(vars("lat", "long", "distance"), funs(as.numeric(as.character(.)))) %>%
  left_join(.,
            data.frame(
              id = kec@data$id,
              kec_code = kec@data$CC_3),
            by = "id")


spcoord = 
  distances %>%
  dplyr::select("lat", "long")

coordinates(spcoord) = ~long + lat

#finding communities in banned area
banned =
  over(spcoord, op) %>%
  mutate(rowid = rownames(.),
         rowid = as.numeric(rowid)) %>%
  mutate(id = case_when(is.na(id) ~ 0,
                        TRUE ~ 1)) %>%
  set_colnames(c("opium_ban", "rowid"))

kec_distances = 
  distances %>%
  rowid_to_column() %>%
  left_join(., banned, by = "rowid") %>%
  mutate(opium_legal = case_when(opium_ban == 0 ~ 1,
                                 opium_ban == 1 ~ 0,
                                 TRUE ~ NA_real_),
         forcing = case_when(opium_legal == 0 ~ (distance - (distance*2)),
                             opium_legal == 1 ~ distance,
                             TRUE ~ NA_real_)) %>%
  dplyr::select(-opium_ban) %>%
  
  #bringing in the mosques
  left_join(.,
            mosques_data %>%
              filter(year_construct < 1810) %>%
              group_by(kec_code = as.character(kec_code)) %>%
              summarise(num_mosques = n()),
            by = "kec_code") %>%
  mutate(num_mosques_nona = case_when(is.na(num_mosques) ~ 0,
                                      TRUE ~ as.numeric(num_mosques)))

            
test_mosques =
  kec_distances %>%
  filter(distance < 15000)

lm(num_mosques_nona ~ opium_legal+ forcing + opium_legal*forcing, data = test_mosques) %>% summary() #test estimate

mosque_plot = loess_plot(x = test_mosques$forcing, W = test_mosques$opium_legal, y = test_mosques$num_mosques_nona, ylab = "Number of Mosques (Pre-1809)", scale = T)

mosque_plot <- mosque_plot + coord_cartesian(ylim = c(0, 2)) # scaling to deal with outliers near threshold
ggsave(plot = mosque_plot, "./_4_outputs/figures/appendix_2e_balance_plot.pdf", height = 3, width = 4)


